www.gusucode.com > 超声波测量以及形成图像 对相关信号进行模拟仿真 > 超声波测量以及形成图像 对相关信号进行模拟仿真/digital holograpy/prog/spsrr.m
function [ varargout ] = spsrr( x0,y0,dx,N,z,varargin )
%SPSRR Single Point Source Recording and Reconstruction
% The complex field formed by a point source is sampled
% on a plane. The sampling interval is dx, N*N is the number
% of sampling points. To ensure the field of the point source
% be sampled adequately, the distance between the point source
% and the sampling plane must be smaller than a minimum value:
% zmin. zmin is decided by :
% zmin=(2*dx*max(x0,y0)+dx^2*N)/lambda
% in which (x0,y0) is the lateral position of the point source.
% the distance between the point source and the sampling plane
% is z, to ensure adequate sampling.,z must be bigger than zmin.
% zmin is an output of the function
%
% Syntax:
% [...,zmin]=spsrr( x0,y0,dx,N,z,'PropertyName','PropertyValue',... );
% Urcs is the reconstructed field on the source plane
% Urcd is the recorded field on the recording plane
%
%--------------------------------------------------------------------------
% PropertyName and PropertyValue:
%
% output {both} | rcs | rcd
% it specifies the style of output
% rcs - the reconstructed field on the object plane is outputted,
% but the recorded field on the sampling plane is not
% outputted.
% [Urcs,dxo,zmin]=spsrr( x0,y0,dx,N,z,'output','rcs' );
% rcd - the recorded field on the sampling plane is outputted,
% but the reconstructed field on the object plane is not
% outputted. and dxo is also not outputted.
% [Urcd,zmin]=spsrr( x0,y0,dx,N,z,'output','rcd' );
% both - both of them and dxo are outputted.
% [Urcs,Urcd,dxo,zmin]=spsrr( x0,y0,dx,N,z,'output','both' );
%
% fa {on} | off
% it decides whether fresnel approximation is applied while obtaining
% the wave field on the recording plane
% on - fresnel approximation applied
% off - fresnel approximation not applied
%
% rm {fresnel} | ASP
% it decides which reconstruction method is used
% fresnel - fresnel transform method
% ASP - Angular Spectrum Propagation method
%
% enlarge {no} | em
% em - if the reconstruction method is 'fresnel', the Urcd is pasted
% on a em*N by em*N all-zero array, if the reconstruction method
% is 'frequency', the frequency array in the 'fdif' function
% (reconstruction function) is pasted on such an all-zero array.
% In this way, details of the PSF can be seen.
% no - not enlarged
%
% inclination {off} | k | rs |
% it decides whether and which inclination factor is considered while
% obtaining the wave field on the recording plane
% k - use Kirchhoff inclination factor
% rs - use Rayleigh-Sommerfeld inclination factor
% off - inclination factor is not considered
%
% fillfactor {off} | on | [ffx,ffy,sampling quality]
% off - the fillfactor of the CCD tends to zero, the sampling of CCD
% is ideal
% on - the fillfactor is taken into consideration in the recording
% simulation, the optical field in each pixel of the CCD is
% sampled on many points, the values of each pixel are evaluated
% by the sum of the values on these points. The sampling
% interval is decided by the average spacial frequency of the
% optical intensity field in each pixel. The fill factor is set
% to be one.
% [ffx,ffy,sampling quality] - specify the fillfactor and the
% sampling quality. If we denote the sampling quality as SQ,
% then the sampling density is SQ times of the default value
% which just satisfies the Nyquist sampling theory
% ------------------------------------------------------------------------
% Urcd is the Point Spread Function of in-line PSDH,
% without considering the integral effect of CCD
% ------------------------------------------------------------------------
error(nargchk(5,17,nargin))
%for n=1:length(varargin)
% if ~ischar(varargin{n})
% error('Property names and values must be characters')
% end
%end
% ----------------------construct the property array----------------------
parray=[0,0,0,0,0,0]; % initialize the property array, the property array
% is used to store the property specified by the user
property=struct('name',{'output','fa','rm',... % construct the property structure
'enlarge','inclination','fillfactor'});
property(1).value={'both','rcs','rcd',''}; % which stores all possible properties
property(2).value={'on','off',''};
property(3).value={'fresnel','ASP',''};
property(4).value={'no',''};
property(5).value={'k','rs','off',''};
property(6).value={'off','on'};
em='1';
for l=1:2:length(varargin)
namefound=0;
valuefound=0;
m=0;
n=0;
while ~namefound && m<length(parray)
m=m+1;
if strcmp(varargin{l},property(m).name)
namefound=1;
end
end
while namefound && ~valuefound && n<length(property(m).value)
n=n+1;
if strcmp(varargin{l+1},property(m).value{n})
valuefound=1;
parray(m)=n;
end
end
%----------------------------------------------
if namefound==1 && m==4 && valuefound==0
parray(4)=3;
em=varargin{l+1};
valuefound=1;
end
%----------------------------------------------
if namefound==0
error('wrong property name')
elseif valuefound==0 && m~=6
error('wrong property value')
end
if namefound==1 && m==6 && valuefound==0
ffp=varargin{l+1};
parray(6)=-1;
end
end
% ------------------------------------------------------------------------
% ------------------check the number of output arguments------------------
% ------------------corresponding to the output property------------------
if parray(1)==0 || parray(1)==1 || parray(1)==4
error(nargchk(4,4,nargout))
elseif parray(1)==2
error(nargchk(3,3,nargout))
elseif parray(1)==3
error(nargchk(2,2,nargout))
end
% ------------------------------------------------------------------------
% ----------calculate the recorded field on the recording plane-----------
load lambda;
zmin=(2*dx*max(x0,y0)+dx^2*N)/lambda;
[vx,vy]=opticimage([N,N],dx);
[X,Y]=meshgrid(vx,vy);
if parray(6)==0 || parray(6)==1
if parray(2)==2
R=sqrt((X-x0).^2+(Y-y0).^2+z.^2);
Urcd=exp(i*2*pi/lambda*R)./R;
else
Urcd=exp(i*2*pi/lambda*z).*exp(i*pi/lambda/z*((X-x0).^2+(Y-y0).^2))./z;
end
if parray(5)==1
ifactor=0.5*(1+z./R);
Urcd=Urcd.*ifactor;
elseif parray(5)==2
ifactor=z./R;
Urcd=Urcd.*ifactor;
end
else
if parray(6)==2
ffp=[1,1,1];
end
snx=dx.*ffp(1).*2.*abs(X-x0)./sqrt((X-x0).^2+z.^2)./lambda;
sny=dx.*ffp(2).*2.*abs(Y-y0)./sqrt((Y-y0).^2+z.^2)./lambda;
sn=ceil(max(snx,sny));
sn=max(sn(:)).*ffp(3);
clear snx sny
I1=zeros(N);
I2=I1;
Io=I1;
cn=1;
if sn>1
[xip,yip]=opticimage([sn,sn],dx*ffp(1)/floor(sn/2)/2,dx*ffp(2)/floor(sn/2)/2);
else
xip=0;
yip=0;
end
[XIP,YIP]=meshgrid(xip,yip);
XIP=repmat(XIP,[1,1,cn*N]);
YIP=repmat(YIP,[1,1,cn*N]);
for n=0:N/cn-1
Xsub=zeros(1,1,cn*N);
Ysub=Xsub;
Xsub(1,1,:)=X(n*cn*N+(1:cn*N));
Ysub(1,1,:)=Y(n*cn*N+(1:cn*N));
Xsub=repmat(Xsub,[sn,sn])+XIP;
Ysub=repmat(Ysub,[sn,sn])+YIP;
R=sqrt((Xsub-x0).^2+(Ysub-y0).^2+z.^2);
Urcdpart=exp(i*2*pi/lambda.*R)./R.*0.5.*(1+z./R);
I1(n*cn*N+(1:cn*N))=sum(sum(abs(1+Urcdpart).^2,1),2)/sn.^2;
I2(n*cn*N+(1:cn*N))=sum(sum(abs(exp(i*pi/2)+Urcdpart).^2,1),2)/sn.^2;
Io(n*cn*N+(1:cn*N))=sum(sum(abs(Urcdpart).^2,1),2)/sn.^2;
end
clear R xip yip XIP YIP Xsub Ysub R Urcdpart
Urcd=[(I1-Io-ones(N))+i*(I2-Io-ones(N))]/2;
clear I1 I2 Io
end
%-----------------------------------------------------------------------
%----------calculate the reconstructed field on the source plane----------
if parray(1)~=3
if parray(3)==2
Urcs=ASP(Urcd,-z,[dx,dx],'pfz',em);
dxo=dx/str2num(em);
else
if em~='1'
[Urcs,dxo]=fresnel(paste(zeros(str2num(em)*N),Urcd),-z,dx);
else
[Urcs,dxo]=fresnel(Urcd,-z,dx);
end
end
end
%-------------------------------------------------------------------------
if parray(1)==0 || parray(1)==1 || parray(1)==4
varargout{1}=Urcs;
varargout{2}=Urcd;
varargout{3}=dxo;
varargout{4}=zmin;
elseif parray(1)==2
varargout{1}=Urcs;
varargout{2}=dxo;
varargout{3}=zmin;
elseif parray(1)==3
varargout{1}=Urcd;
varargout{2}=zmin;
end